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1.  Introduction 


Over  the  past  decade,  the  use  of  improvised  explosive  devices  (IEDs)  by  enemy 
forces  has  increased  exposure  among  US  combat  vehicles  and  Soldiers  to  under¬ 
body  blast  (UBB)  events.  During  such  an  event,  explosive  detonated  underground 
causes  a  rapid  expansion  of  soil  and  gases  under  a  vehicle.  Momentum  is  imparted 
to  the  vehicle  hull  by  a  combination  of  the  expanding  gases  and  the  surrounding 
soil.  Energy  is  transmitted  through  the  hull  to  the  internal  vehicular  floor  that  within 
milliseconds  produces  a  deformation  of  the  floor  plate.  These  rapid  deformations 
of  the  floor  plate  transmit  significant  loads  to  the  Soldier’s  lower  extremities  and 
induce  accelerative  injury.1,2 

Medical  injury  data  reveals  that  the  lower  extremities  are  one  of  the  most  com¬ 
monly  injured  body  regions  among  mounted  Soldiers  in  theater.3^6  The  current  un¬ 
derstanding  of  these  lower  extremity  injuries  is  largely  based  on  research  within 
the  car  crash  industry.  However,  UBB  events  exert  forces  on  occupants  of  ground 
combat  vehicles  that  are  of  higher  magnitude,  in  different  directions,  and  over  a 
shorter  duration  than  forces  in  civilian  car  accidents,1,6  making  these  injuries  and 
their  mechanisms  unique  to  the  military  population. 

During  Operation  Enduring  Freedom  and  Operation  Iraqi  Freedom  more  than  16,000 
Soldiers  were  unable  to  return  to  duty  due  to  their  injuries.7  Recent  studies  have 
shown  that  up  to  50%  of  Soldiers  wounded  in  action  have  sustained  a  lower  ex¬ 
tremity  injury.7^9  UBB-induced  injuries  to  the  lower  extremities  are  frequently  as¬ 
sociated  with  fractures  to  multiple  bones  and  may  also  involve  damage  to  the  sur¬ 
rounding  soft  tissues,  nerves,  and  vascular  network.  It  has  been  shown  that  for  pa¬ 
tients  that  have  sustained  multiple  injuries,  those  with  foot  and  ankle  injuries  have 
a  greater  disability  score  than  those  without  foot  and  ankle  injuries.  Furthermore, 
those  who  sustained  an  injury  to  the  calcaneus  are  associated  with  the  highest  com¬ 
plication  rate.8  Advances  in  technology  have  both  increased  the  opportunities  for 
surgeons  to  salvage  limbs,  as  well  as  helped  in  the  development  of  better  prosthe- 
ses  for  amputees.  Patients  who  decide  to  salvage  their  lower  extremities  often  re¬ 
quire  a  longer  period  of  rehabilitation  compared  to  those  who  decide  to  amputate 
their  leg.10  To  the  wounded  warrior,  a  return  to  high  levels  of  physical  activity  is  of¬ 
ten  a  priority  since  it  will  affect  their  ability  to  return  to  full  military  duty.  It  is  not 
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uncommon  for  patients  with  salvaged  limbs  who  have  dealt  with  chronic  infection, 
constant  pain,  or  poor  functional  outcome  to  opt  for  elective  amputation.710  Under¬ 
standing  the  mechanism  that  causes  one  injury  over  another  is  therefore  important 
for  protection  design  concepts. 

The  process  of  insult  to  injury  must  be  well  understood  to  effectively  design  protec¬ 
tion  systems  that  mitigate  lower  extremity  injuries  of  mounted  Soldiers  induced  by 
IEDs.  Physical  models,  including  human  cadavers,  and  frangible  and  mechanical 
test  devices,  are  used  to  develop  a  better  understanding  of  injury  mechanisms  in  the 
lower  extremities.31 1-15  Human  cadaver  models  are  often  accepted  as  the  best  rep¬ 
resentation  of  the  living  human  leg  due  to  their  high  level  of  biofidelity.  The  draw¬ 
back  of  human  cadaver  models,  however,  is  that  there  are  significant  anatomical 
variations  between  cadavers,  including  anthropometry  and  tissue  properties.  Fur¬ 
thermore,  the  age  and  health  status  of  the  cadaver  at  time  of  death  may  not  be  rep¬ 
resentative  of  the  Soldier  population  and  it  is  largely  unknown  how  well  the  proper¬ 
ties  of  cadaveric  tissue  represents  that  of  living  tissue.  Frangible  leg  models  are  de¬ 
signed  to  be  representative  of  an  “average”  Soldier  and  eliminate  the  variation  in  an¬ 
thropometry  and  tissue  properties  between  specimens.  However,  the  material  prop¬ 
erties  are  homogeneous  and  not  necessarily  accurate  simulants  of  biological  tissue. 
In  addition,  frangible  models  generally  can  only  be  used  once  and  therefore  are  not 
considered  cost  effective  for  extensive  testing.  In  contrast,  mechanical  test  devices, 
such  as  anthropomorphic  test  devices  (ATDs)  can  survive  multiple  tests  since  they 
are  constructed  out  of  rigid  materials.  However,  the  rigidness  of  the  model  leads  to 
a  different  mechanical  response  than  the  more  compliant  human  at  high  rates  of  ac¬ 
celeration.  Thus,  the  stiff  mechanical  response  of  ATDs  is  not  always  suitable  for 
studying  injury  mechanisms  of  humans  in  UBB  events. 

Numerical  models  are  the  alternative  to  the  physical  models  discussed  above.  The 
strength  of  a  numerical  model  is  that  it  can  easily  be  changed  to  explore  the  role 
and  importance  of  a  physical  effect.  These  models  can  simulate  a  wide  range  of 
UBB  conditions  on  the  lower  extremities  and  capture  the  location  of  where  injury 
initiates  and  how  it  propagates.  Thus  they  can  evaluate  the  effectiveness  of  protec¬ 
tion  concepts  such  as  energy-absorbing  structures.  However,  as  with  any  numeri¬ 
cal  tool,  numerical  models  necessarily  make  assumptions  on  geometries,  boundary 
conditions,  constitutive  response,  and  material  parameters. 
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Numerically  modeling  the  human  is  a  particularly  challenging  effort  requiring  care¬ 
ful  attention  to  detail.  Numerical  models  used  to  study  UBB  events  in  the  literature 
have  largely  been  adapted  from  mature  models  in  the  car  crash  injury,  however, 
this  approach  may  not  be  appropriate  due  to  the  unique  loading  of  UBB  events. 
The  large  variability  in  human  anthropometry  and  experimental  uncertainty  in  tis¬ 
sue  characterization  is  largely  underrepresented  in  the  simulation  literature.  In  fact, 
only  a  small  subset  of  typical  material  models  and  their  parameters  are  used  to  nu¬ 
merically  model  the  human  body  in  UBB  events  and  the  anthropometry  of  the  hu¬ 
man  model  is  held  fixed. 

Biological  variability,  including  tissue  properties  and  the  anthropometric  differ¬ 
ences  between  Soldiers,  may  play  a  critical  role  in  injury  mechanisms.  Numerical 
understanding  of  injury  mechanisms  is  simplified  by  starting  with  a  simpler  prob¬ 
lem  and  building  physical  insight.  Here  we  present  a  baseline  study  that  only  evalu¬ 
ates  the  role  of  material  parameters  on  injury  mechanisms.  The  report  is  organized 
as  follows.  Section  2  discusses  the  details  of  the  model  and  the  boundary  conditions 
of  the  simulations.  Section  3  presents  the  results  of  our  study  including  an  overview 
(Section  3.1),  a  detailed  look  at  2  illustrative  simulations  (Section  3.2),  the  trends 
of  predicted  injuries  (Section  3.3),  the  effects  of  specific  parameters  on  potential  in¬ 
jury  occurrence  (Section  3.4),  and  the  location  of  predicted  failures  (Section  3.5). 
Finally,  Section  4  discusses  our  observations  of  the  results  and  considers  the  future 
direction  of  the  model  and  its  applications. 

2.  Simulation  Details  and  Injury  Prediction 

This  section  describes  the  anatomical  detail  captured  in  our  numerical  model,  the 
applied  loading  condition,  and  our  choices  in  constitutive  models.  Our  method  for 
determining  whether  a  bone  would  have  likely  failed  due  to  the  applied  loading  is 
also  discussed. 
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2.1  Anatomical  Details 


Figure  1  provides  an  overview  of  the  simulation  geometry,  anatomical  details,  and 
the  loading  condition.  Initially  developed  in  Kraft  et  al., 16  the  lower  leg  model  was 
modified  to  produce  the  finite  element  model  used  in  this  work  (see  Appendix  A  for 
additional  details).  Figure  la  shows  the  flesh  layer  (yellow),  cortical  bone  (gray), 
trabecular  bone  (orange,  see  inset)  and  the  steel  plate  (black)  used  in  the  simula¬ 
tions.  Thirty  bones  are  explicitly  represented,  including  the  femur,  patella,  tibia, 
fibula,  and  the  26  bones  of  the  foot.  Both  cortical  and  trabecular  bone  are  modeled 
in  the  femur,  patella,  tibia,  fibula,  talus,  and  calcaneus.  The  bones  of  the  mid-  and 
forefoot  are  not  segmented  into  different  bone  types  and  are  assigned  cortical  bone 
properties.  The  soft  tissue  and  skin,  also  referred  to  as  flesh,  is  a  homogenization 
of  skin,  muscles,  tendons,  connective  tissue,  blood  vessels,  etc.  using  an  approach 
similar  to  that  proposed  by  Cheung  et  al.17  Note  that  there  are  more  than  100  mus¬ 
cles,  tendons,  and  ligaments  in  the  lower  leg  that  are  simplified  in  taking  this  as¬ 
sumption.  A  flesh  layer  fills  in  the  gaps  between  the  bones  and  thus  details  of  the 
joint  structures  such  as  cartilage  and  synovial  fluid  are  not  explicitly  represented. 
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Time  (ms) 


Fig.  1  Panel  a  provides  an  overview  of  the  geometry  considered  in  our  simulations.  The  inset 
in  panel  a  (dashed  red  box)  shows  a  cross-sectional  view  of  the  talus  and  calcaneus  where  the 
cortical  shell  (gray)  and  trabecular  bone  (orange)  are  shown  against  the  flesh  (yellow).  Panel  b 
shows  a  bone  map,  a  2-dimensional  (2-D)  representation  of  the  bones  of  the  lower  leg  excluding 
the  patella  and  femur.  Labels  are  provided  for  all  bones.  Note  that  this  presentation  does  not 
maintain  all  of  the  spatial  connectivity,  e.g.,  the  calcaneus  is  not  in  direct  contact  with  the  tibia. 
Panel  c  is  a  trace  of  the  velocity  versus  time  that  is  applied  to  the  steel  plate  in  the  positive  Z- 
direction. 
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The  anatomical  geometry  for  the  model  was  obtained  from  Zygote  Media  Group, 
Inc.  (American  Fork,  UT),  a  company  that  creates  anatomical  geometry  from  CT 
scans  of  human  specimens  who  anthropometrically  fit  within  the  50th  percentile. 
The  femur  and  patella  were  positioned  so  that  the  lower  leg  is  in  a  seated  position 
(90°  angles  at  both  the  knee  and  ankle)  and  initialized  so  that  the  foot  rests  on  a 
steel  plate. 

2.2  Loading  and  Boundary  Conditions 

An  idealized  velocity  history  (Fig.  lc)  is  imposed  on  the  steel  plate,  which  is  treated 
as  rigid.  This  loading  condition  is  similar  to  those  seen  in  blast  events,  but  is  not 
meant  to  reproduce  a  specific  threat,  merely  to  capture  gross  features.  Special  care  is 
taken  to  ensure  the  applied  velocity  and  its  derivative  are  continuous,  thus  avoiding 
spurious  noise.18  The  mathematical  formulation  of  the  smooth  curve  is  described 
in  Appendix  B .  Initially,  the  plate  is  stationary  and  then  over  the  course  of  2  ms  is 
accelerated  to  a  constant  velocity  of  5.5  m/s  before  accelerating  back  to  zero.  After 
10  ms,  the  plate  comes  to  rest  with  a  displacement  of  3  cm.  In  the  first  2  ms,  the 
average  acceleration  of  the  plate  is  2,750  m/s2,  thus  acceleration  due  to  gravity  is 
ignored.  The  femur  is  not  fixed  in  any  direction. 

2.3  Material  Details 

Extensive  testing  on  bone  and  soft  tissues  reveal  rate-dependencies  and  nonlineari¬ 
ties  in  their  constitutive  responses.  Accurately  predicting  fracture  patterns  of  bone 
and  subsequent  injury  requires  a  rigorous  understanding  of  the  constitutive  and  fail¬ 
ure  response  of  the  material  as  well  as  adequate  representations  of  interfaces  and 
contact  surfaces — a  resolution  that  is  beyond  the  current  state  of  the  model.  To  pro¬ 
vide  a  baseline  to  compare  against,  we  take  the  rate  dependence  and  nonlinearity  of 
the  materials  as  secondary  in  importance  and  assume  that  the  flesh,  cortical  bone, 
and  trabecular  bone  behave  as  isotropic  linear  elastic  materials.  As  we  will  discuss 
in  Section  2.4,  we  use  state  variables  as  correlative  and  predictive  of  injury  but  do 
not  attempt  to  model  the  actual  failure  process.  We  assume  that  injury  occurs  from 
accumulated  strain,  and  thus  we  treat  any  rate-dependent  failure  to  be  of  secondary 
importance  in  the  problem. 
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Table  1  is  a  compilation  of  constitutive  models  and  their  parameters  used  to  repre¬ 
sent  tissues  of  the  lower  leg  found  within  the  literature.  The  first  column  gives  the 
anatomic  component  and  the  second  column  lists  the  types  of  constitutive  models 
used.  The  third  column  provides  the  range  of  material  properties  specific  to  a  linear 
elastic  constitutive  model  that  we  analyzed  in  this  study.  In  cases  where  the  con¬ 
stitutive  model  was  not  linear  elastic,  the  range  is  determined  by  either  taking  the 
small  strain  elastic  modulus  or  an  equilibrium  modulus.  The  last  column  provides 
the  references  for  the  materials.  After  surveying  the  literature  for  a  range  of  mate¬ 
rial  properties  used  to  model  tissues  in  the  lower  leg,  we  selected  6  parameters  that 
we  varied  between  relatively  low  or  high  values  that  fell  within  the  ranges  reported 
in  the  literature.  The  anatomical  components  and  the  material  model  parameters 
that  we  vary  are  reported  in  Table  2.  For  example,  in  a  given  simulation  the  Young’s 
modulus  of  cortical  bone  was  either  15  or  19  GPa.  The  density,  Young’s  modulus 
and  Poisson’s  ratio  of  trabecular  bone,  and  the  bulk  modulus  and  Poisson’s  ratio  of 
flesh  are  also  varied.  Since  the  flesh  is  treated  as  linear  elastic,  the  reporting  of  its 
bulk  modulus  can  unambiguously  be  replaced  with  a  Young’s  modulus.  In  all  of  the 
simulations  the  density  and  Poisson’s  ratio  of  cortical  bone,  as  well  as  the  density 
of  flesh,  remained  fixed.  To  handle  the  flesh  material  parameter  values  numerically, 
we  utilize  a  nodal-based  tetrahedron,19  see  Appendix  A. 
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Table  1  Compilation  of  various  constitutive  models  and  parameters  found  in  computational 
models  of  the  lower  extremities.  Failure  stress  is  given  in  terms  of  a  von  Mises  stress  criterion. 


Anatomic 

Component 

Material  Model 

Material  Properties 

References 

Trabecular  bone 

Linear  elastic. 
Elastic-plastic, 
Transversely  isotropic 

E  =  150-900  MPa 
v  =  0.3-0.45 
p  =  650-1,100  kg/m3 

1,20-23 

Cortical  bone 

Lineal-  elastic. 
Elastic-plastic, 
Transversely  isotropic 

E  =  15-19  GPa 

v  =  0.3-0.35 

p  =  1,680-2,000  kg/m3 
Failure  Stress  =  124-175 

MPa 

1 ,20-24 

Flesh 

Lineal'  elastic. 
Hyper-elastic, 

Lineal'  viscoelastic 

K  =  2.5-300  MPa 

=  0.42-0.49 
p  =  1,000-1,300  kg/m3 

1,20,21,23 

Note:  Not  all  parameters  are  listed  for  each  constitutive  model.  Furthermore,  not  all  lower  extrem¬ 
ity  components  are  listed.  E  is  the  Young’s  modulus,  v  is  Poisson’s  ratio,  p  is  density,  and  K  is  the 
bulk  modulus. 


Table  2  Linear  elastic  material  parameters.  Three  material  parameters  were  kept  constant 
and  6  material  parameters  were  varied  through  low  or  high  values  to  produce  64  combinations. 

Anatomic 

Material  Model 

Low  Value 

High  Value 

Fixed  Value 

Component 

Parameter 

Cortical  bone 

P 

1,850  kg/m3 

E 

15  GPa 

19  GPa 

V 

0.3 

Trabecular  bone 

P 

650  kg/m3 

1,100  kg/m3 

E 

150  MPa 

900  MPa 

V 

0.3 

0.45 

Flesh 

P 

1,000  kg/m3 

K 

20  MPa 

300  MPa 

V 

0.42 

0.48 

Note:  E  denotes  Young’s  modulus,  v  denotes  Poisson’s  ratio,  p  is  density,  and  K  is  the  bulk  modulus. 


2.4  Correlating  State  Variables  to  Injury  Prediction 

Injury  can  occur  over  multiple  length  scales  and  involve  cellular  processes  as  well 
as  the  mechanical  disruption  of  tissue.  The  numerical  simulations  solve  for  the  me¬ 
chanics  of  the  problem,  but  correlating  a  state  variable  or  variables  with  injury  is  an 
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ongoing  field  of  research.  Typically,  principal  stress  or  strain,  peak  stress  or  strain, 
pressure,  and  von  Mises  stress  are  used  as  the  criteria  to  trigger  failure  in  finite  el¬ 
ement  models  through  element  deletion,  degradation  of  the  elastic  modulus,  or  co¬ 
hesive  failure  in  elements.  However,  once  failure  occurs,  these  numerical  simula¬ 
tions  typically  diverge  from  what  is  seen  experimentally,  as  fracture  is  a  difficult 
phenomenon  to  model.  Accurately  predicting  the  particular  injury  of  bone  fracture 
requires  a  detailed  understanding  of  material  properties  as  well  as  the  microstruc¬ 
ture  of  the  tissue.  Since  the  spatial  resolution  of  our  model  is  insufficient  to  predict 
precise  fracture  patterns,  our  focus  is  whether  certain  bones  would  likely  have  frac¬ 
tured  during  an  event.  These  simulations  establish  trends  of  how  particular  combi¬ 
nations  of  material  parameters  affect  the  observed  stresses. 

We  assume  that  a  mechanical  state  variable  such  as  pressure  or  von  Mises  stress 
can  be  used  as  an  injury  predictor,  i.e.,  a  correlative  measure  for  a  condition  that 
would  ultimately  result  in  injury.  We  are  using  the  convention  that  normal  stresses 
are  negative  in  compression  and  define  pressure  as  negative  one-third  the  trace  of 
the  Cauchy  stress.  We  restrict  ourselves  to  stresses  observed  in  the  cortical  shell 
bone.  The  constitutive  model  we  use  to  represent  the  cortical  tissue  is  isotropic 
linear  elastic  and  implies  that  strain  and  stress  can  be  interchanged  uniquely.  The 
specific  choice  of  an  injury  predictor  is  an  open  research  endeavor.  As  one  might 
imagine,  the  reduction  of  the  time  history  of  a  symmetric  rank-3  tensor  to  an  in¬ 
jury  predictor  is  extremely  open  ended  and  material  specific.  Without  experimental 
measurements  to  conclusively  reveal  a  predictor,  we  use  von  Mises  stress  and  pres¬ 
sure  since  they  are  both  reduced  representations  of  the  complete  stress  state  and  are 
scalar  valued.  The  von  Mises  stress  is  a  typical  indicator  of  failure  in  other  mate¬ 
rials  and  the  ranges  in  Table  1  correspond  to  the  elastic-plastic  failure  values  used 
for  cortical  bone.  Pressure,  on  the  other  hand,  is  not  commonly  used  as  an  indicator 
of  injury  in  simulations.  However,  it  also  represents  a  directionally  free  and  poten¬ 
tially  physically  meaningful  choice. 

Other  common  choices  for  injury  predictors  include  “axial  stress”,  principal  stresses, 
and  the  maximum  of  the  principal  stresses.  The  complicated  geometry,  especially  at 
interfaces,  makes  axial  stress  not  universal  to  the  entire  bone  geometry.  Similarly, 
the  ambiguity  in  using  the  principal  stresses  arises  from  a  complicated  geometry 
and  stress  state,  i.e.,  there  is  little-to-no  meaning  of  the  eigenvectors  from  element 
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to  element  or  from  time  to  time  making  their  values  difficult  to  interpret.  The  no¬ 
tion  of  the  maximum  of  the  eigenvalues  of  ever  changing  directions  also  carries  lit¬ 
tle  practical  meaning.  If  a  microstructural  detail  was  known,  however,  and  included 
within  the  model,  e.g.,  osteon  direction,  the  stress  oriented  with  this  microstructural 
detail  might  be  an  extremely  useful  injury  predictor. 

The  approach  we  take  to  predicting  injury  is  as  follows.  The  completed  simulation 
is  postprocessed  to  determine  whether  any  elements  within  the  cortical  shell  exceed 
an  injury  predictor  threshold.  The  2  injury  predictors  we  consider  are  von  Mises 
stress  and  magnitude  of  pressure.  We  select  a  wide  range  of  threshold  values  so  that 
we  can  assess  any  sensitivity  of  predicted  injury  on  the  specific  choice  of  threshold. 
For  the  von  Mises  stress  injury  predictor  we  consider  the  range  of  failure  reported 
in  Table  1 .  Since  it  is  not  typically  used  in  the  literature,  the  range  of  interest  was 
widened  for  the  pressure  predictor.  The  interpretation  of  the  data  needs  to  be  care¬ 
fully  considered  since  the  transmission  of  stress  through  bone  will  change  between 
fractured  and  unfractured  cases. 

Postprocessing  was  performed  using  GROPE  and  ALGEBRA  (SEACAS  Toolkit, 
Sandia  National  Laboratory),  and  MATLAB.  Using  GROPE  and  ALGEBRA,  el¬ 
ements  can  be  probed  using  text  commands  following  the  completion  of  a  simu¬ 
lation.  Additionally  calculations  can  be  made  quite  rapidly,  including  picking  the 
maximum  value  over  time.  These  tools  were  used  to  determine  the  pressure  and 
von  Mises  stress  over  time  (for  each  microsecond)  for  all  elements. 

3.  Results 

Six  material  parameters  of  the  cortical  bone,  trabecular  bone,  and  flesh  were  var¬ 
ied  between  low  and  high  values,  producing  64  combinations  of  material  proper¬ 
ties.  The  specific  values  of  material  properties  are  reported  in  Table  3.  Each  of  these 
simulations  can  be  thought  of  as  a  proxy  for  a  Soldier  with  a  slightly  different  tissue 
characterization,  all  within  the  accepted  range  of  values  shown  in  Table  1 .  The  re¬ 
sults  are  organized  as  follows.  Section  3.1  provides  an  overview  of  the  simulations. 
Section  3.2,  considers  2  simulations  in  detail.  Section  3.3,  comments  on  the  trends 
of  the  injury  predictor  found  across  the  data  sets.  Section  3.4  reports  on  the  occur¬ 
rence  of  failure  based  on  the  material  parameters  themselves.  Finally,  we  report  on 
the  localization  of  stress  in  Section  3.5. 
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Table  3  Material  properties  of  64  simulations  produced  by  varying  6  material  parameters 


# 

Ecort 

Arab 

Etrab 

t'trab 

KFlesh 

^Flesh 

# 

Ecort 

Arab 

Etrab 

Arab 

KFlesh 

^Flesh 

(GPa)  (kg/m3)  (GPa) 

(MPa) 

(GPa)  (kg/m3)  (GPa) 

(MPa) 

1 

15 

650 

0.15 

0.3 

20 

0.42 

33 

19 

650 

0.15 

0.3 

20 

0.42 

2  / 

15 

650 

0.15 

0.3 

20 

0.48 

34 

19 

650 

0.15 

0.3 

20 

0.48 

3 

15 

650 

0.15 

0.3 

300 

0.42 

35 

19 

650 

0.15 

0.3 

300 

0.42 

4 

15 

650 

0.15 

0.3 

300 

0.48 

36 

19 

650 

0.15 

0.3 

300 

0.48 

5 

15 

650 

0.15 

0.45 

20 

0.42 

37 

19 

650 

0.15 

0.45 

20 

0.42 

6 

15 

650 

0.15 

0.45 

20 

0.48 

38 

19 

650 

0.15 

0.45 

20 

0.48 

7 

15 

650 

0.15 

0.45 

300 

0.42 

39 

19 

650 

0.15 

0.45 

300 

0.42 

8 

15 

650 

0.15 

0.45 

300 

0.48 

40 

19 

650 

0.15 

0.45 

300 

0.48 

9 

15 

1100 

0.15 

0.3 

20 

0.42 

41 

19 

1100 

0.15 

0.3 

20 

0.42 

10 

15 

1100 

0.15 

0.3 

20 

0.48 

42 

19 

1100 

0.15 

0.3 

20 

0.48 

11 

15 

1100 

0.15 

0.3 

300 

0.42 

43 

19 

1100 

0.15 

0.3 

300 

0.42 

12 

15 

1100 

0.15 

0.3 

300 

0.48 

44 

19 

1100 

0.15 

0.3 

300 

0.48 

13 

15 

1100 

0.15 

0.45 

20 

0.42 

45 

19 

1100 

0.15 

0.45 

20 

0.42 

14 

15 

1100 

0.15 

0.45 

20 

0.48 

46 

19 

1100 

0.15 

0.45 

20 

0.48 

15 

15 

1100 

0.15 

0.45 

300 

0.42 

47 

19 

1100 

0.15 

0.45 

300 

0.42 

16 

15 

1100 

0.15 

0.45 

300 

0.48 

48 

/ 

19 

1100 

0.15 

0.45 

300 

0.48 

17 

15 

650 

0.90 

0.3 

20 

0.42 

49 

/ 

19 

650 

0.90 

0.3 

20 

0.42 

18 

15 

650 

0.90 

0.3 

20 

0.48 

50 

/ 

19 

650 

0.90 

0.3 

20 

0.48 

19 

15 

650 

0.90 

0.3 

300 

0.42 

51 

/ 

19 

650 

0.90 

0.3 

300 

0.42 

20 

15 

650 

0.90 

0.3 

300 

0.48 

52 

19 

650 

0.90 

0.3 

300 

0.48 

21 

15 

650 

0.90 

0.45 

20 

0.42 

53 

19 

650 

0.90 

0.45 

20 

0.42 

22 

15 

650 

0.90 

0.45 

20 

0.48 

54 

/ 

19 

650 

0.90 

0.45 

20 

0.48 

23 

15 

650 

0.90 

0.45 

300 

0.42 

55 

19 

650 

0.90 

0.45 

300 

0.42 

24 

15 

650 

0.90 

0.45 

300 

0.48 

56 

/ 

19 

650 

0.90 

0.45 

300 

0.48 

25 

15 

1100 

0.90 

0.3 

20 

0.42 

57 

19 

1100 

0.90 

0.3 

20 

0.42 

26 

15 

1100 

0.90 

0.3 

20 

0.48 

58 

19 

1100 

0.90 

0.3 

20 

0.48 

27 

15 

1100 

0.90 

0.3 

300 

0.42 

59 

19 

1100 

0.90 

0.3 

300 

0.42 

28 

15 

1100 

0.90 

0.3 

300 

0.48 

60 

19 

1100 

0.90 

0.3 

300 

0.48 

29 

15 

1100 

0.90 

0.45 

20 

0.42 

61 

19 

1100 

0.90 

0.45 

20 

0.42 

30 

15 

1100 

0.90 

0.45 

20 

0.48 

62 

/ 

19 

1100 

0.90 

0.45 

20 

0.48 

31 

15 

1100 

0.90 

0.45 

300 

0.42 

63 

19 

1100 

0.90 

0.45 

300 

0.42 

32 

15 

1100 

0.90 

0.45 

300 

0.48 

64 

/ 

19 

1100 

0.90 

0.45 

300 

0.48 

Note:  A  check  mark  indicates  that  the  simulation  terminated  early,  i.e.,  before  3.5  ms.  Blue  rows 
indicate  the  2  illustrative  simulations  discussed  in  Section  3.2.  The  4  green  rows  isolate  a  specific 
subset  where  flesh  properties  are  varied  and  are  discussed  in  Section  3.3.  The  7  orange  rows  are  a 
subset  where  predicted  injuries  followed  a  different  trend  discussed  in  Section  3.3. 
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3.1  Summary  of  Simulations 


In  total,  64  simulations  were  executed  for  all  the  combinations  of  material  param¬ 
eters  discussed  in  Section  2.  These  simulations  are  enumerated  as  rows  in  Table  3. 
Each  column  entry  in  Table  3  gives  the  value  of  a  material  parameter  used  in  that 
particular  simulation.  As  described  in  Section  2.2,  every  simulation  experiences  the 
same  loading  condition.  Cauchy  stress  tensor  data  was  output  for  all  elements  every 
0.1  [is  (simulation  time).  Each  simulation  used  128  cores  and  was  allowed  a  wall¬ 
time  of  48  hours.  Typically,  this  walltime  was  sufficiently  long  enough  for  a  simula¬ 
tion  to  reach  17  ms  (simulation  time)  and  output  400  GB  of  raw  data  per  simulation. 
Regular  hard  drive  write  commands  slowed  the  performance  of  the  solver.  Of  the  64 
simulations,  55  of  them  reached  or  exceeded  3.5  ms  of  simulation  time  before  ei¬ 
ther  terminating  from  a  negative  volume  or  running  out  of  walltime  due  to  a  plum¬ 
meting  time  step.  Simulations  that  failed  before  3.5  ms  are  denoted  in  Table  3  by  a 
checkmark  under  the  column  labeled  ET.  Thirteen  rows  in  Table  3  are  highlighted 
blue,  green  or  orange.  The  blue  rows  are  2  illustrative  simulations  (Simulations  1 
and  3)  that  have  been  separated  out  and  studied  in  more  detail  in  Section  3.2.  The  4 
green  rows  (Simulations  9  through  12)  and  the  7  orange  rows  (Simulations  18,  22, 
26,  30,  50,  54,  and  58)  are  discussed  in  Section  3.3. 

The  55  simulations  that  reached  or  exceeded  3.5  ms  represent  roughly  86%  of  the 
total  simulations.  By  3.5  ms,  the  plate  is  displaced  6  mm  and  reaches  a  velocity  of 
4.6  m/s,  which  is  84%  of  the  5.5  m/s  maximal  value  of  the  velocity  loading  seen  in 
Fig.  lc.  Only  26  simulations  (40%)  reached  or  exceeded  5  ms  of  simulation  time. 
The  plate  in  these  simulations  reaches  maximal  velocity  and  has  the  largest  dis¬ 
placements  (>14  mm).  Of  the  9  simulations  that  terminated  before  3  ms  of  simula¬ 
tion  time  (denoted  by  a  check  mark  in  Table  3),  7  of  them  had  a  similar  set  of  mate¬ 
rial  properties.  These  include  a  cortical  Young’s  modulus  of  19  GPa  and  a  trabecular 
Young’s  modulus  of  900  MPa.  Five  of  the  9  simulations  that  terminated  early  had  a 
cortical  Young’s  modulus  of  19  GPa,  a  trabecular  Young’s  modulus  of  900  MPa,  and 
a  trabecular  density  of  650  kg/m3.  Seven  of  the  9  simulations  that  terminated  early 
had  the  higher  Poisson’s  ratio  for  the  soft  tissue  of  0.48.  It  is  unclear  whether  these 
parameter  combinations  are  significant  for  numerical  stability.  However,  none  of  the 
simulations  with  the  low  soft  tissue  bulk  modulus  exceeded  5  ms  of  simulation  time. 
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3.2  Two  Illustrative  Cases 


This  section  considers  2  simulations  in  detail.  The  only  difference  in  material  pa¬ 
rameters  is  the  bulk  modulus  of  the  flesh.  The  material  properties  corresponding  to 
these  simulations  are  highlighted  in  blue  in  Table  3,  where  Simulation  1  had  the  low 
soft  tissue  bulk  modulus  value  (20  MPa)  and  Simulation  3  had  the  high  soft  tissue 
bulk  modulus  value  (300  MPa).  The  remaining  material  properties  in  these  simula¬ 
tions  were  assigned  low  values.  We  first  present  time  history  results  and  then  con¬ 
sider  the  injury  prediction  capabilities  of  these  simulations,  i.e.,  the  sensitivity  and 
susceptibility  of  specific  bones  to  failure. 

Figures  2  and  3,  respectively,  show  the  time  histories  of  the  pressure  and  von  Mises 
stress  localization  in  the  deformed  configuration.  The  2  perspectives  were  chosen 
to  highlight  some  of  the  structurally  important  regions  of  the  calcaneus  and  talus. 
Each  row  corresponds  to  a  different  time  indicated  at  the  far  left  side  of  the  figure. 
The  results  from  Simulation  1  (low  bulk  modulus)  are  presented  in  the  left  column 
and  the  results  from  Simulation  3  (high  bulk  modulus)  are  presented  in  the  right 
column.  Note  that  the  maximum  values  on  the  color  scales  are  quite  different  for 
Figs.  2  and  3,  i.e.,  different  magnitudes  were  seen  when  comparing  von  Mises  stress 
and  pressure. 
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Simulation  3  70  MPa 


-70  MPa 
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Simulation  1 


Fig.  2  Pressure  localization  throughout  time  for  Simulation  1  and  Simulation  3.  Two  different 
views  of  the  model  are  presented  for  each  simulation.  The  time  of  the  simulation  is  indicated 
at  the  far  left  side  of  the  figure  and  applies  to  both  simulations.  The  *  in  the  figure  indicates 
locations  that  exceed  the  color  scale. 
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Fig.  3  von  Mises  stress  localization  throughout  time  for  Simulation  1  and  Simulation  3.  Two 
different  views  of  the  model  are  presented  for  each  simulation.  The  time  of  the  simulation  is 
indicated  at  the  far  left  side  of  the  figure  and  applies  to  both  simulations. 
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The  time  history  results  for  pressure  are  qualitatively  similar  between  simulations. 
As  time  progresses,  the  magnitude  of  pressure  increases  and  the  stress  states  are 
similar  between  simulations,  i.e.,  compressive  regions  and  tensile  regions  appear 
in  similar  anatomical  locations.  As  one  looks  at  the  false  color  throughout  time, 
the  color  patterns  only  increase  in  intensity  as  opposed  to  changing  in  color.  In 
both  simulations,  compressive  pressures  are  seen  at  the  interface  between  the  cal¬ 
caneus  and  talus,  and  the  pressure  is  largely  tensile  on  the  inferior  surface  of  the 
calcaneus  (on  what  would  be  considered  the  “under  side”  of  the  calcaneus)  nearest 
the  cuboid  and  cuneiforms.  In  addition  to  the  similarities  between  the  2  calcanei, 
there  are  similarities  in  the  pressure  distribution  for  the  fibula  and  talus  as  well.  The 
fibula  experiences  compressive  pressure  just  above  the  ankle  joint  (lateral  malleo¬ 
lus).  The  talus  experiences  compressive  pressure  where  it  meets  the  calcaneus  and 
cuneiforms.  There  are  similarities  in  the  stress  localization  when  considering  the 
von  Mises  stress  as  well.  Figure  3  shows  that  the  calcaneus  again  experiences  large 
stresses  on  its  inferior  surface.  Furthermore,  there  is  some  localization  of  the  shear 
stresses  in  the  fibula. 

Quantitatively,  however,  the  2  simulations  are  quite  different.  At  0.1  ms,  the  pres¬ 
sures  in  Simulation  3  can  be  seen  as  a  light  pink  (Fig.  2).  In  Simulation  1,  these 
pressures  are  not  readily  visible  until  the  2.0  ms  panel.  Four  milliseconds  into  Sim¬ 
ulation  3,  the  portion  of  the  calcaneus  nearest  the  steel  plate  (calcaneal  tuberosity) 
experiences  compressive  pressures  that  exceed  the  color  scale,  whereas  in  Simula¬ 
tion  1  they  are  significantly  smaller.  Similar  observations  can  be  made  when  com¬ 
paring  the  color  intensity  values  in  Fig.  3. 

Single  element  time  histories  of  the  pressure  and  von  Mises  stress  from  the  cortical 
shell  of  6  bones  are  presented  in  Fig.  4.  Each  color  corresponds  to  a  different  bone 
and  the  relative  location  of  each  element  is  indicated  with  a  circle  in  panel  a.  An 
open  circle  indicates  that  the  element  is  on  the  side  of  the  bone  not  visible  from 
the  current  view  point.  Pressure  results  from  Simulations  1  and  3  are  plotted  in 
panels  b  and  c,  respectively,  and  von  Mises  stress  results  are  shown  in  panels  d 
and  e.  The  inset  in  panels  b-e  show  the  early  time  (less  than  1  ms)  behavior  of  those 
same  points.  The  vertical  dashed  lines  in  panels  b  and  d  indicate  the  time  at  which 
Simulation  1  terminated. 
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Fig.  4  Time  histories  for  the  pressure  and  von  Mises  stress  from  selected  elements  in  the 
cortical  layer  of  various  bones.  Each  color  corresponds  to  an  element  in  a  particular  bone. 
Pressure  results  from  Simulations  1  and  3  are  plotted  in  panels  b  and  c  and  von  Mises  stress  in 
panels  d  and  e.  The  inset  in  each  figure  shows  the  early  time  behavior. 
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Simulation  3  provides  a  clear  overview  of  the  result  of  the  loading  on  the  lower  leg. 
Panel  c  presents  the  pressure  histories  throughout  time.  For  early  times,  i.e.,  times 
less  than  1  ms,  the  elements  oscillate  between  compressive  and  tensile  states  (see 
inset).  This  is  most  notable  in  the  calcaneus  (red  trace).  During  this  initial  tran¬ 
sient  state,  the  amplitude  of  the  oscillations  are  small  with  values  less  than  0.5%  of 
the  peak  values  observed  in  the  plot.  For  times  greater  than  1  ms,  the  transient  os¬ 
cillatory  behavior  gives  way  to  a  smooth  accumulation  of  pressure  followed  by  a 
smooth  reduction  of  pressure.  The  largest  peak  pressure  occurs  in  the  fibula  (black 
trace,  panel  c),  which  reaches  a  maximum  value  of  approximately  30  MPa  around 
4.2  ms.  This  is  followed  by  the  tibia  (blue  trace)  with  a  pressure  of  approximately 
20  MPa  around  4.3  ms  and  the  calcaneus  (red  trace)  with  a  pressure  of  15  MPa  at 
4.1  ms.  During  the  times  presented,  the  element  in  the  talus  is  the  only  one  to  expe¬ 
rience  a  tensile  (negative)  pressure.  Qualitatively  similar  observations  can  be  made 
for  the  von  Mises  stresses  in  Simulation  3  (panel  e).  During  the  initial  transient 
period,  there  are  high-frequency  oscillations  in  von  Mises  stress  of  the  calcaneus. 
The  elements  from  the  other  bones,  however,  exhibit  fewer  oscillations.  The  initial 
transient  period  gives  way  to  a  smooth  accumulation  of  stress  that  rises  to  a  peak 
value  and  then  falls.  The  von  Mises  stress  peaks  largest  in  the  calcaneus  at  100  MPa 
around  4.1  ms,  followed  by  the  fibula  at  90  MPa  at  approximately  the  same  time. 

The  early  time  behavior  for  the  pressure  (panel  b,  inset)  in  Simulation  1  is  also  os¬ 
cillatory  in  nature  and  small  in  amplitude.  However,  comparing  the  calcaneus  traces 
between  the  2  simulations  reveals  higher  frequency  content  present  in  the  signal  for 
Simulation  3,  the  high  case  for  the  bulk  modulus  of  the  flesh.  The  differences  in 
these  2  oscillations  is  likely  due  to  ringing  occurring  in  the  soft  tissue  of  the  heel 
pad  between  the  calcaneus  and  the  steel  plate. 

After  the  initial  transient  period,  both  the  pressure  and  von  Mises  stress  in  Simu¬ 
lation  1  increase  smoothly  up  to  the  point  at  which  the  simulation  terminates.  Al¬ 
though  impossible  to  verify,  based  on  the  nature  of  the  loading,  one  might  assume 
that  the  qualitative  features  of  Simulation  1  would  be  similar  to  Simulation  3.  We 
can,  however,  compare  the  peak  values  between  panels  b  and  c  or  stress  in  panels  d 
and  e.  In  Simulation  3,  the  pressure  in  the  fibula  (black  trace,  panel  c)  reaches  a 
maximum  value  of  approximately  30  MPa  around  4.1  ms.  Around  this  same  time, 
the  fibula  in  Simulation  1  (black  trace,  panel  b),  has  a  pressure  of  about  15  MPa.  A 
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similar  observation  can  be  made  by  comparing  the  fibula  probe  in  panels  d  and  e, 
however,  in  this  case  Simulation  1  reaches  approximately  70%  of  the  peak  value 
experienced  in  Simulation  3. 

At  4.1  ms  it  is  clear  that  the  distribution  of  stress  is  different  between  the  2  simula¬ 
tions.  For  instance,  compare  the  order  of  colored  traces  from  largest  to  smallest  in 
panels  d  and  e.  In  Simulation  1  (panel  d)  the  order  from  largest  stress  to  smallest 
is  the  probe  in  the  fibula,  calcaneus,  talus,  tibia,  cuboid,  and  femur.  In  Simulation  3 
(panel  d)  the  order  is  calcaneus,  fibula,  tibia,  femur,  talus,  and  cuboid.  These  are  re¬ 
sults  from  single  element  probes  but  whose  values  strongly  depend  on  the  value  of 
the  bulk  modulus  of  the  flesh. 


After  the  initial  transient  period,  the  time  histories  presented  in  Fig.  4  show  a  clear 
maximum  value  for  both  pressure  and  von  Mises  stress.  Thus,  we  reduce  the  data 
further  by  considering  peak  values  of  pressure  and  von  Mises  stress,  which  we 
present  in  a  2-D  schematic  of  the  lower  leg  or  bone  map  (Figs.  5  and  6).  These  fig¬ 
ures  allow  a  quick  assessment  of  a  particular  bone’s  sensitivity  to  reaching  injuri¬ 
ous  pressures  or  von  Mises  stress  levels  during  a  simulation  and  thus  can  be  used  to 
assess  trends  in  predicted  injury. 
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Fig.  5  Bone  maps  using  the  pressure  magnitude  criteria  for  Simulations  1  and  3  (highlighted 
in  blue  in  Table  3) 
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Fig.  6  Bone  maps  using  the  von  Mises  criteria  for  Simulations  1  and  3  (highlighted  in  blue  in 
Table  3).  The  blue  dashed  box  represents  the  failure  range  as  considered  in  the  current  study. 


Figures  5  and  6  compare  the  pressure  or  von  Mises  injury  predictor  for  Simulation  1 
and  Simulation  3,  respectively.  Each  column  of  the  figure  corresponds  to  a  particular 
threshold  level.  A  bone  that  is  colored  red  indicates  that  at  least  one  element  in  the 
cortical  layer  of  that  bone  exceeded  the  threshold,  whereas  a  gray  bone  signifies  that 
no  elements  of  that  bone  exceeded  the  threshold.  If  a  bone  has  been  highlighted  red, 
the  number  of  elements  that  exceeded  the  threshold  is  noted  on  or  near  the  bone. 
The  numbers  provide  some  insight  into  the  significance  of  the  color  coding,  i.e.,  if 
only  a  single  element  exceeds  the  threshold  it  is  unlikely  to  be  a  strong  indicator. 
The  bones  of  the  forefoot  are  grouped  as  one  bone  for  this  analysis.  A  naming  key 
for  the  bones  can  be  found  in  Fig.  lb. 

Figure  5  presents  the  pressure  criteria  used  for  predicting  injuries.  Similar  to  the 
observations  made  in  Fig.  4,  the  stresses  in  Simulation  3  were  typically  larger  than 
those  in  Simulation  1.  The  first  column  shows  that  the  tibia,  calcaneus,  and  talus  in 
Simulation  1  experience  pressure  larger  in  magnitude  than  20  MPa.  The  second  and 
third  columns  show  that  only  the  talus  meets  the  noted  threshold.  The  number  of 
elements  that  exceed  the  thresholds  in  these  cases  are  most  significant  for  the  talus 
and  calcaneus.  In  contrast,  the  first  column  for  Simulation  3  shows  that  each  bone 
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has  at  least  one  element  that  exceeds  20  MPa  of  pressure.  The  numbers  of  elements 
in  this  case  are  significant  for  all  of  the  bones  except  the  forefoot.  The  2  bones  that 
experience  the  largest  pressures  in  Simulation  3  are  the  calcaneus,  which  exceeds 
100  MPa,  and  the  talus,  which  exceeds  60  MPa. 

Figure  6  presents  the  von  Mises  injury  predictor  for  Simulation  1  and  3.  In  con¬ 
trast  to  Fig.  5,  the  magnitude  of  stress  in  Fig.  6  is  much  larger.  The  range  of  failure 
thresholds  used  in  the  literature  and  summarized  in  Table  1  is  indicated  in  Fig.  6  by 
the  dashed  blue  rectangle.  The  von  Mises  injury  predictor  for  Simulation  1  shows 
that  the  talus  exceeds  a  von  Mises  stress  of  260  MPa,  but  that  no  other  bones  ex¬ 
ceed  80  MPa.  Simulation  3,  in  contrast,  shows  that  both  the  calcaneus  and  talus  ex¬ 
ceed  260  MPa,  the  largest  threshold  exceeded  by  the  fibula  was  140  MPa,  followed 
by  the  tibia,  which  reaches  100  MPa.  Recall  that  the  only  difference  between  Simu¬ 
lation  1  and  Simulation  3  was  the  value  of  the  bulk  modulus  for  the  flesh.  For  these 
2  simulations,  the  higher  soft  tissue  bulk  modulus  led  to  larger  stresses  in  the  corti¬ 
cal  shells  of  the  bones  within  the  observation  time  frame. 

3.3  Trends  in  Predictive  Injury 

Here  we  present  the  trends  observed  after  reducing  the  data  sets  to  their  respective 
bone  maps.  Eleven  specific  simulations  are  also  discussed,  which  are  highlighted  in 
either  green  or  orange  in  Table  3. 

After  surveying  both  the  von  Mises  and  the  pressure  injury  predictor  bone  maps 
for  all  of  the  simulations,  we  made  the  following  observations.  Using  either  injury 
predictor,  the  most  commonly  injured  bones  were  the  talus  and  calcaneus.  For  a 
fixed  failure  threshold,  very  few  injuries  were  predicted  in  the  forefoot  compared 
to  the  other  bones.  It  was  common  to  see  both  the  talus  and  the  calcaneus  share 
the  same  injury  prediction  sensitivity,  but  there  were  some  interesting  deviations 
from  this  trend.  In  some  cases,  the  calcaneus  and  talus  were  so  susceptible  to  injury 
that  even  the  largest  von  Mises  threshold  of  260  MPa  was  exceeded.  This  is  seen 
in  Simulations  3,  7,  11,  15,  35,  39,  43,  and  47,  which  all  ran  with  the  higher  flesh 
bulk  modulus  and  the  lower  trabecular  bone  Young’s  modulus  value.  The  only  cases 
where  the  calcaneus  was  more  susceptible  to  injury  than  the  talus  were  cases  where 
we  analyzed  the  pressure.  This  occurred  for  Simulations  3,  7,  8,  11,  15,  16,  35,  36, 
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39,  40,  43,  44,  47,  and  48,  which  is  essentially  the  same  set  of  simulations  as  the 
group  where  the  calcaneus  and  talus  exceeded  the  260  MPa  von  Mises  threshold. 


In  Simulations  9-12,  only  the  Poisson’s  ratio  or  bulk  modulus  of  flesh  was  varied 
(green  rows,  Table  3).  In  these  simulations,  the  properties  for  cortical  and  trabecular 
bone  were  assigned  the  same  low  or  high  value.  The  von  Mises  injury  predictor  and 


pressure  injury  predictor  bone  maps  are  presented  in  Figs.  7  and  8,  respectively.  For 
each  figure,  a  row  corresponds  to  a  different  simulation. 
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Fig.  7  Bone  maps  using  the  von  Mises  threshold  criteria  for  Simulations  9-12  (highlighted 
green  in  Table  3).  The  blue  dashed  box  represents  the  failure  range  as  considered  in  the  current 
study. 
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Fig.  8  Bone  maps  using  the  pressure  threshold  criteria  for  Simulations  9-12  (highlighted  green 
in  Table  3) 


Figure  7  considers  the  sensitivity  to  bone  failure  for  the  4  simulations  as  it  depends 
on  the  von  Mises  stress  injury  predictor.  The  only  difference  between  Simulations 
9  and  10,  and  Simulations  11  and  12  is  the  value  of  Poisson’s  ratio.  Based  on  the 
number  of  elements  that  have  reached  a  threshold  criteria  in  the  highlighted  bones 
(red),  we  predict  that  Poisson’s  ratio  had  little  effect  on  the  amount  of  damage  we 
see  in  the  bones.  The  parameter  that  has  a  larger  role  in  the  number  of  injured 
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bones  is  the  bulk  modulus  of  flesh.  This  observation  is  made  by  comparing  Simu¬ 
lations  9  and  10,  which  were  run  with  the  lower  bulk  modulus  value,  against  Sim¬ 
ulations  1 1  and  12,  which  were  run  with  the  higher  value.  Consider,  for  example,  a 
single  von  Mises  stress  threshold  of  140  MPa.  In  Simulation  9,  there  are  6  elements 
in  the  talus  that  meet  this  criteria,  and  for  Simulation  10  there  is  only  one.  In  con¬ 
trast,  the  talus  in  Simulations  11  and  12  had  18  and  20  elements,  respectively,  that 
met  the  criteria.  Simulation  11  also  showed  that  3  elements  in  the  fibula  and  113 
elements  in  the  calcaneus  reached  the  injury  threshold,  and  Simulation  12  showed 
that  5  elements  in  the  calcaneus  reached  this  same  injury  threshold. 

By  considering  a  single  bone  and  varying  the  threshold,  we  also  get  an  intuitive 
picture  of  the  sensitivity  of  injury  prediction  for  that  particular  bone.  The  talus  in 
Simulations  9,  11,  and  12  predicts  injury  regardless  of  the  threshold,  thus  making 
it  virtually  insensitive  to  our  injury  predictor.  This  was  a  common  trend  observed 
seen  for  the  talus.  In  contrast,  predicting  injury  in  the  calcaneus  is  sensitive  both  to 
the  threshold  of  the  injury  predictor  as  well  as  the  material  properties  used  in  the 
simulations. 

Figure  8  presents  the  predicted  injuries  for  Simulations  9  through  12  using  the 
pressure  injury  criteria.  Save  for  the  lowest  injury  threshold  of  20  MPa  in  Simu¬ 
lation  11,  the  forefoot  region  does  not  show  any  elements  above  threshold  for  in¬ 
jury  prediction  and  is  therefore  virtually  insensitive  to  the  pressure  injury  criteria  in 
these  4  simulations.  The  other  bones,  however,  do  show  some  sensitivity,  although 
for  much  lower  thresholds  than  the  von  Mises  criteria.  The  most  interesting  feature 
that  comes  from  analyzing  the  pressure  injury  predictor  is  seen  by  comparing  Sim¬ 
ulations  9  and  10  with  Simulations  1 1  and  12.  Recall  that  using  the  von  Mises  stress 
injury  predictor,  the  talus  was  virtually  insensitive  to  the  threshold  value  and  typ¬ 
ically  predicted  failure.  Using  the  pressure  injury  predictor  we  begin  to  see  cases 
where  the  talus  is  no  longer  overly  sensitive  to  predicting  injury.  By  comparing  the 
talus  and  the  calcaneus  in  Simulations  9  and  10,  one  can  see  that  the  talus  experi¬ 
ences  larger  pressures  than  the  calcaneus.  However,  this  trend  is  reversed  for  Sim¬ 
ulations  11  and  12  where  the  calcaneus  experiences  larger  pressures  than  the  talus. 
This  feature  is  important  since  it  implies  that  a  pressure  injury  criteria  failure  thresh¬ 
old  can  be  picked  such  that  the  calcaneus  will  predict  failure  without  the  talus  pre¬ 
dicting  failure  as  well. 
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We  conclude  the  discussion  of  the  bone  map  results  by  considering  Simulations  18, 
22,  26,  30,  50,  54,  and  58,  the  orange  rows  in  Table  3.  These  are  cases  where  the 
talus  was  extremely  susceptible  to  pressure  predicted  injury,  but  the  calcaneus  was 
not,  i.e.,  the  talus  was  always  above  threshold  and  the  calcaneus  was  rarely  predicted 
as  injured.  Of  these  7  simulations,  5  reached  simulation  times  of  3.5  ms  but  did  not 
greatly  exceed  3.5  ms,  and  the  remaining  2  only  reached  3.1  ms.  All  7  simulations 
had  the  same  trabecular  Young’s  modulus  of  900  MPa,  soft  tissue  properties  of  bulk 
modulus  20  MPa,  and  Poisson’s  ratio  of  0.48. 

3.4  Distribution  of  Predicted  Injury 

Here  we  survey  the  occurrence  of  predicted  injuries  specific  to  individual  bones.  As 
reported  in  Section  3.1,  the  talus  and  calcaneus  were  common  locations  of  injury. 
The  histograms  in  Figs.  9  and  10  report  the  general  effect  of  a  material  parameter  on 
predicted  injury.  Failure  occurrence  is  reported  for  the  calcaneus,  talus,  tibia,  fibula 
and  forefoot  based  on  the  value  used  for  that  material  parameter  (Figs.  9,  10,  and 
Appendix  C).  Failure  is  said  to  occur  during  a  simulation  for  a  particular  bone  if  at 
least  one  element  in  the  cortical  shell  exceeds  the  injury  predictor  threshold.  Both 
von  Mises  stress  and  pressure  injury  predictors  were  considered.  These  histograms 
were  created  for  each  of  the  6  material  parameters,  and  the  2  failure  criteria. 


Calcaneus  Talus  Tibia  Fibula  Forefoot  Key 


Fig.  9  Failure  occurrence  in  each  bone  varying  the  bulk  modulus  of  flesh  using  a  pressure 
injury  predictor.  Failure  occurrences  are  tallied  at  various  bones  for  a  range  of  thresholds. 
Each  panel  is  specific  to  an  individual  bone  and  tallies  the  number  of  simulations  for  which  the 
absolute  value  of  the  pressure  within  any  element  in  the  bone  exceeded  the  damage  threshold. 
The  bar  colors  indicate  which  damage  threshold  the  tally  corresponds  to.  Low  and  high  values 
for  the  material  parameter  as  well  as  the  specific  damage  thresholds  used  are  indicated  on  the 
right  under  the  key. 
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Fig.  10  Failure  occurrence  in  each  bone  varying  the  bulk  modulus  of  flesh  using  a  von  Mises 
stress  injury  predictor.  Failure  occurrences  are  tallied  at  various  bones  for  a  range  of  thresh¬ 
olds.  Each  panel  is  specific  to  an  individual  bone  and  tallies  the  number  of  simulations  for 
which  the  von  Mises  stress  within  any  element  in  the  bone  exceeded  the  damage  threshold. 
The  bar  colors  indicate  which  damage  threshold  the  tally  corresponds  to.  Low  and  high  val¬ 
ues  for  the  material  parameter  as  well  as  the  specific  damage  thresholds  used  are  indicated  on 
the  right  under  the  key. 


Figure  9  divides  the  simulations  by  their  soft  tissue  bulk  modulus  and  considering 
the  results  from  using  the  pressure  injury  predictor.  Histograms  are  shown  for  the 
calcaneus,  talus,  tibia,  fibula,  and  forefoot.  Failure  occurrence  refers  to  a  maximum 
of  32  simulations,  in  either  the  low  or  high  category.  Seven  threshold  values  were 
considered  to  investigate  the  sensitivity  of  failure.  These  7  values  include  and  extend 
beyond  the  range  listed  in  Table  1 .  The  value  of  the  threshold  is  indicated  by  the 
color  of  the  bar  in  the  key.  The  key  also  contains  the  number  of  simulations  that 
did  not  reach  3.5  ms  and  therefore  were  omitted  from  the  histogram.  Consider,  for 
example,  the  panel  dedicated  to  the  calcaneus  in  Fig.  9.  The  tall,  dark  blue  bar  that 
reaches  the  number  20  indicates  that  20  out  of  the  32  simulations  that  were  run  with 
a  high  flesh  bulk  modulus  had  at  least  one  element  exceed  the  pressure  damage 
threshold  of  80  MPa.  Similarly,  Fig.  10  looks  at  the  effect  of  the  flesh  bulk  modulus 
on  von  Mises  failure.  Figures  C-l  through  C-10  in  Appendix  C  explore  the  same 
question  for  the  other  5  material  parameters  that  were  varied  in  the  study. 

The  material  property  that  had  the  strongest  influence  on  occurrences  of  predicted 
injury  was  the  bulk  modulus  of  the  flesh.  This  effect  is  most  clearly  demonstrated  in 
Fig.  10.  By  considering  the  panel  dedicated  to  the  calcaneus,  one  can  see  that,  when 
using  the  high  bulk  modulus  value  of  300  MPa,  a  significant  number  of  simulations 
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predicted  calcaneus  failures  across  the  different  damage  thresholds  (e.g.,  28  simu¬ 
lations  out  of  32  for  a  threshold  value  of  120  MPa).  There  is  also  a  high  occurrence 
of  predicted  injury  in  the  tibia  and  fibula  in  the  high  bulk  modulus  cases.  In  con¬ 
trast,  the  low  bulk  modulus  cases  do  not  predict  any  injuries  in  the  calcaneus,  tibia, 
or  fibula.  The  Poisson’s  ratio  of  the  soft  tissue  and  the  Young’s  modulus  of  cortical 
and  trabecular  bone  exhibit  a  minor  effect  on  predicted  injury.  There  was  no  clear 
trend  of  predicted  injury  on  the  density  or  Poisson’s  ratio  of  trabecular  bone.  Note 
that  none  of  the  low  cases  for  the  soft  tissue  bulk  modulus  exceed  a  simulation  time 
of  5  ms. 

3.5  Stress  Localization 

Sections  3.1  and  3.4  reported  that  injuries  to  the  calcaneus  and  talus  were  common 
amongst  the  simulations.  While  the  histograms  presented  in  Section  3.4  show  the 
number  of  simulations  with  these  injuries,  the  specific  element  location  is  not  re¬ 
ported.  Here  we  consider  the  localization  of  the  largest  stresses  of  the  problem. 

Each  panel  in  Fig.  1 1  presents  the  calcaneus  from  1  of  6  different  orthogonal  views. 
The  calcaneus  has  been  colored  dark  gray  with  additional  colors  superimposed  on 
top  of  it.  The  color  indicates  whether  that  specific  element  in  at  least  1  of  the  64  sim¬ 
ulations  exceeded  a  specific  threshold,  thus  blocks  of  uniform  color  indicate  multi¬ 
ple  elements  from  multiple  simulations  all  meeting  the  same  criteria.  In  each  panel, 
the  color  red  indicates  pressures  larger  than  80  MPa,  blue  represents  pressures  less 
than  -80  MPa,  and  green  signifies  von  Mises  stress  exceeding  150  MPa.  These  3 
colors  are  treated  as  additive  so  that  the  appearance  of  yellow  in  the  panel  means 
that  both  the  green  and  red  conditions  have  been  meet.  Cyan  illustrates  that  the  blue 
and  green  conditions  have  been  met,  and  black  indicates  all  3  conditions  have  been 
met.  An  additive  of  red  and  blue,  i.e.,  both  pressure  extremes  without  a  von  Mises 
extreme  was  not  observed.  The  large  red  patch  at  the  calcaneal  tuberosity,  most  vis¬ 
ible  in  panels  d,  e,  and  /  reveals  a  patch  of  elements  experiencing  large  compres¬ 
sive  pressures.  This  observation  was  also  made  specifically  for  Simulations  1  and  3, 
when  discussing  Fig.  2.  This  patch  of  elements  are  on  the  portion  of  the  calcaneus 
closest  to  the  steel  plate.  The  large  green  patch,  most  visible  in  panels  a  and  e  corre¬ 
sponds  to  a  region  of  large  von  Mises  stress  concentration  along  the  medial  side  of 
the  calcaneus.  The  last  notable  collection  of  elements  are  multicolored  and  can  be 
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seen  in  panel  e,  a  view  through  the  phalanges  to  the  calcaneus.  Both  large  von  Mises 
stresses  and  tensile  pressures  are  observed  on  the  inferior  surface  of  the  calcaneus, 
the  point  on  the  calcaneus  nearest  the  talus. 


Fig.  11  Stress  concentrations  of  the  calcaneus.  Red  indicates  pressures  larger  than  80  MPa 
were  experienced  at  that  location  for  a  least  1  of  the  64  simulations.  Blue  indicates  pressures 
lower  than  -80  MPa.  Green  indicates  von  Mises  stresses  larger  than  150  MPa.  The  colors  are 
additive,  therefore  yellow  means  that  both  the  green  and  red  conditions  have  been  met,  cyan 
indicates  that  the  blue  and  green  conditions  have  been  met,  and  black  is  representative  of  all  3 
conditions.  Each  panel  shows  a  different  view  of  the  foot. 


Figure  12  is  the  analogous  figure  for  showing  the  stress  concentrations  for  the  talus. 
Panels  b,  d,  and /  show  multiple  patches  of  von  Mises  stress  concentrations  (green) 
and  combinations  of  other  extremes  (blue,  yellow,  and  black)  in  and  around  the 
joint  of  the  2  bones.  Panels  c  and  e  reveal  additional  patches  of  elements  where  the 
von  Mises  stresses  are  large  and  the  pressures  are  compressive.  These  regions  are 
located  near  the  calcaneus  but  also  where  the  talus  meets  the  tibia. 
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Fig.  12  Stress  concentrations  of  the  talus.  Red  indicates  pressures  larger  than  80  MPa  were 
experienced  at  that  location  for  a  least  1  of  the  64  simulations.  Blue  indicates  pressures  lower 
than  -80  M  Pa.  Green  indicates  von  Mises  stresses  larger  than  150  MPa.  The  colors  are  ad¬ 
ditive  therefore  yellow  means  that  both  the  green  and  red  conditions  have  been  met,  cyan  in¬ 
dicates  that  the  blue  and  green  conditions  have  been  met,  and  black  is  representative  of  all  3 
conditions.  Each  panel  shows  a  different  view  of  the  foot. 
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4.  Discussion 


We  conducted  64  simulations  on  a  lower  leg  finite  element  model  to  understand  the 
sensitivity  of  predicted  injury  on  the  choices  of  material  parameters.  The  same  ac¬ 
celerative  loading  condition  was  applied  to  each  of  the  64  simulations.  This  work  is 
one  of  many  studies  to  be  conducted  on  the  lower  leg  and  was  designed  to  serve  as  a 
baseline  for  future  comparisons.  Here  we  discuss  relevant  issues  that  have  not  been 
addressed  previously  in  the  report.  For  example,  several  significant  unanswered  re¬ 
search  questions  are  as  follows:  What  physical  effects  are  most  important  for  pre¬ 
dicting  injury?  Does  injury  occur  on  a  strain  and  time  scale  where  nonlinearity  and 
rate-dependencies  are  large  effects,  or  can  these  effects  be  neglected?  Furthermore, 
have  we  assumed  a  sufficient  level  of  detail  in  our  model? 

4.1  Biofidelity  and  Injury  Predictive  Capability 

Throughout  the  simulation  literature,  there  is  a  large  range  of  tissue  material  prop¬ 
erties.  Here  we  considered  both  low  and  high  values  reported  in  the  literature  for 
6  different  material  parameters.  Using  combinations  of  these  parameters,  we  re¬ 
produced  situations  that  would  be  potentially  injurious  for  bones  of  the  lower  leg. 
While  we  do  not  actually  model  failure,  there  are  some  clear  cases  where  the  con¬ 
centration  of  stress  would  result  in  element  deletion  or  trigger  some  other  numerical 
method  for  failure.  In  Sections  3.1  and  3.4,  we  reported  that  the  soft  tissue  plays  a 
role  in  the  types  of  injuries  predicted.  Two  questions  that  still  need  to  be  addressed: 
Why  does  the  soft  tissue  property  directly  affect  the  stress  concentrations  we  see  in 
the  cortical  shell  of  a  neighboring  bone?  Have  we  represented  the  interfaces  ade¬ 
quately  in  our  numerical  model? 

From  our  results,  we  can  conclude  that  injury  prediction  by  itself  is  not  a  meaningful 
metric  for  validation  purposes.  Figures  5  and  6  illustrate  that  by  only  changing 
one  material  parameter  we  could  significantly  affect  the  specific  bones  that  were 
injured  and  their  severity.  Thus,  certain  combinations  of  material  parameters  could 
selectively  produce  a  different  set  of  injuries.  The  conundrum  lies  in  that  all  of  the 
material  parameters  we  simulated  were  within  the  accepted  range  of  values,  thus 
each  one  represented  a  valid  choice  in  the  model  and  each  one  predicted  different 
injuries. 
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Predicted  calcaneus  and  talus  injuries  were  a  core  focus  of  our  results.  In  Sec¬ 
tion  3.1,  we  reported  that  the  talus  was  highly  susceptible  to  both  pressure  and 
von  Mises  stress  concentrations.  This  was  additionally  supported  through  the  his- 
torgrams  (Figs.  9  and  10,  and  supplemental  figures  in  Appendix  C  where  the  num¬ 
ber  of  talus  failures  far  exceeded  calcaneus  failures  using  a  von  Mises  failure  crite¬ 
ria,  and  are  on  par  with  calcaneus  failures  using  the  pressure  criteria.  Stress  local¬ 
izations  were  identified  in  Figs.  11  and  12  that  pointed  to  the  particular  interface 
where  injuries  typically  would  occur.  Similar  fractures  of  the  talus  were  also  ob¬ 
served  in  Kraft  et  al., 16  whose  model,  despite  sharing  the  same  foundation  as  the 
current  model,  contained  significant  differences  as  well.  For  example,  the  Kraft  et 
al.  model  included  ligaments,  rate-dependent  constitutive  models,  fracture  criteria 
for  cortical  and  trabecular  bone,  and  considered  an  entirely  different  loading  func¬ 
tion.  Despite  these  large  differences,  Kraft  et  al.  also  observed  talus  fractures,  even 
in  cases  where  the  calcaneus  fractured  quite  extensively. 

However,  in  comparison  with  other  lower  extremity  injuries,  fractures  of  the  talus 
are  relatively  uncommon  and  are  typically  the  result  of  high-energy  trauma.  Talus 
fractures  constitute  less  than  1%  of  all  reported  fractures  and  only  3%  to  6%  of 
all  foot  fractures.25-26  Henderson  et  al.  conducted  a  study  in  which  18  postmortem 
human  subject  (PMHS)  lower  limbs  were  subjected  to  a  high  rate  impact,  similar 
to  loads  observed  in  military  vehicles.  In  all  of  their  tests  they  only  saw  3  talus 
fractures.3  The  accepted  injury  mechanism  of  talus  fractures  is  a  hyperdorsiflexion 
force,  in  which  the  toes  are  forced  toward  the  shin  beyond  their  normal  range.  This 
causes  the  posterior  ligaments  of  the  subtalar  joint  (meeting  point  of  the  talus  and 
calcaneus)  to  rupture,  the  neck  of  the  talus  to  impact  against  the  leading  anterior 
edge  of  the  distal  tibia,  and  a  fracture  line  to  develop.25,26 

Due  to  the  unrealistic,  high  number  of  talus  failures  we  saw  in  our  results,  and  the 
localization  of  the  largest  stresses  occurring  between  the  talus  and  the  calcaneus, 
we  believe  that  the  talus  injuries  represent  a  false-positive  in  our  injury  predictive 
capability.  It  is  necessary  to  revisit  and  improve  the  level  of  detail  in  our  model, 
specifically  in  the  complex  joint  structures  of  the  ankle.  Since  calcaneus  injuries 
were  closely  related  to  the  soft  tissue  properties  and  both  talus  and  calcaneus  in¬ 
juries  existed  in  the  original  model  by  Kraft  et  al.,  we  suspect  the  oversight  in  how 
the  joint  is  modeled  is  common  to  both  models.  As  mentioned  in  Section  2,  we 
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represented  all  the  muscles,  tendons,  ligaments,  connective  tissue,  and  skin  as  one 
material.  This  soft  tissue  is  in  a  tied  contact  with  the  bones.  Early  in  the  research, 
shared  nodes  were  used  between  the  soft  tissue  and  the  cortical  bone.  These  simu¬ 
lations  terminated  extremely  quickly  due  to  negative  volumes,  potentially  indicat¬ 
ing  a  mesh  sensitivity  that  was  only  reduced  by  the  definition  of  contact.  In  the  hu¬ 
man  body,  there  are  undoubtedly  both  tied  and  slip  surfaces  that  exist  between  bone 
and  cartilage,  bone  and  ligaments/tendons,  and  bone  and  muscle.  From  a  mechanics 
standpoint,  the  interfaces  between  soft  tissues  and  bone  can  be  quite  challenging  to 
model.  True  joint  structures  transition  from  bone  to  cartilage  to  a  fluid  layer  before 
meeting  the  neighboring  cartilage  layer.  It  is  possible  that  the  current  contact  defi¬ 
nition  between  the  cortical  bone  and  the  soft  tissue  is  overly  restrictive  and  part  of 
the  source  of  large  shear  stresses  in  what  one  might  think  is  a  largely  compression- 
dominated  problem.  We  suspect  the  homogenization  and  perhaps  (overly)  restric¬ 
tive  contact  definitions  are  likely  candidates  to  revisit  and  improve  in  our  model  be¬ 
fore  further  use. 

4.2  Physical  Considerations 

In  this  section,  we  provide  some  estimates  using  isotropic  linear  elastic  theory  that 
we  later  use  to  frame  further  discussion  of  our  results.  Using  the  small  strain  theory 
of  linear  elasticity,  one  can  estimate  the  time  it  takes  for  an  elastic  wave  to  propagate 
from  the  foot  to  the  femoral  head.  The  largest  acoustic  wave  speed  originates  from 
the  cortical  bone  and  is  approximately  3,700  m/s  (for  the  theory  of  linear  elasticity 
see,  for  example,  Landau  and  Lifshitz27).  The  distance  from  the  calcaneus  to  the 
femoral  head  is  approximately  0.9  m,  yielding  an  estimated  travel  time  of  0.24  ms. 
One  can  then  loosely  think  of  0.24  ms  as  a  unit  of  time  that  it  takes  the  wave  front 
to  traverse  the  model  (understanding  of  course  that  there  is  no  continuous  pathway 
of  cortical  bone  from  the  calcaneus  to  the  femur).  The  longest  dimension  of  the 
calcaneus  is  approximately  10  times  smaller  (0.07  m),  and  the  time  it  takes  a  wave 
to  travel  this  distance  is  0.02  ms.  Finally,  the  typical  thickness  of  cortical  bone 
can  vary  but  is  between  millimeters  and  centimeters,  thus  placing  the  travel  times 
on  the  order  of  multiple  microseconds.  In  3.5  ms,  an  elastic  wave  in  cortical  bone 
can  travel  a  distance  of  13  m,  which  is  about  14  lengths  of  the  leg,  190  calcaneus 
lengths,  and  up  to  thousands  of  cortical  thicknesses.  These  estimates  suggest  that 
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at  the  millisecond  time-scale,  elastic  waves  have  ample  time  to  traverse  the  cortical 
shell  thickness  and  the  small  bones  of  the  feet  numerous  times. 

Wave  propagation  within  the  cortical  shell  provides  insight  into  only  part  of  the 
problem.  Additionally,  the  tissues  of  the  body  can  be  thought  of  in  terms  of  layered 
media.  Brekhovskikh’s  book28  treats  a  number  of  problems  in  linear  elasticity  for 
layered  media.  Using  simple  calculations,  one  can  obtain  qualitative  details  regard¬ 
ing  reflection  and  transmission  coefficients  for  layered  tissues  of  the  body.  First, 
consider  the  much- simplified  problem  of  acoustic  waves  in  2  semi-infinite  planes. 
At  normal  incidence,  and  taking  the  same  material  parameter  ranges  used  in  our 
simulations,  we  can  estimate  a  range  in  the  magnitude  of  reflection  coefficients  be¬ 
tween  the  flesh  and  cortical  bone  to  be  0.82-0.96.  Similarly  for  trabecular  and  cor¬ 
tical  bone,  we  find  a  range  of  magnitudes  of  the  reflection  coefficients  to  be  0.60- 
0.86.  These  numbers  provide  some  additional  insight  into  the  effects  the  material 
parameters  can  have  on  transmission  and  reflection  of  elastic  waves.  Next,  consider 
a  wave  that  propagates  from  a  semi-infinite  volume  of  steel  followed  by  a  finite 
thickness  of  flesh  and  then  into  a  semi-infinite  volume  of  cortical  bone.  Calculat¬ 
ing  the  frequency  dependence  of  the  transmission  coefficient  at  normal  incidence 
into  the  bone  illustrates  the  importance  of  understanding  the  layered  tissues  of  our 
problem,  since  the  thickness  of  the  flesh  layer  heavily  influences  the  transmission 
coefficient.  Assuming  for  this  example  a  flesh  thickness  of  6  mm,  we  make  the  fol¬ 
lowing  observations  about  the  transmission  coefficient.  For  frequencies  from  c o  =0 
to  lu  =0.1  MHz,  the  transmission  coefficient  rapidly  falls  from  a  value  of  1  to  al¬ 
most  0.  For  frequencies  larger  than  0. 1  MHz,  the  periodicity  of  the  solution  for  the 
transmission  coefficient  makes  the  magnitude  of  the  transmission  coefficient  cy¬ 
cle  between  complete  transmission  and  complete  suppression.  Thus,  the  thin  flesh 
layer  can  influence  the  transmitted  wave  in  a  similar  way  to  that  of  a  thin  film  in 
optics.  Further  considerations  should  be  made  to  understand  other  multilayer  prob¬ 
lems,  e.g.,  steel  to  flesh  to  cortical  bone  to  trabecular  bone  or  steel  to  rubber  to  flesh 
to  cortical  bone  to  trabecular  bone,  since  it  is  these  thin  layers  that  can  largely  in¬ 
fluence  the  transmitted  wave  and  thus  influence  the  predicted  injuries. 

The  mechanism  for  injury  needs  to  be  better  understood  to  properly  assign  the  ap¬ 
propriate  level  of  detail  to  the  finite  element  model.  Regardless  of  the  mechanism, 
fracture  of  the  bone  tissue  will  occur  within  the  micro  structure  and  the  accumu- 
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lation  of  which  will  have  effects  on  the  tissue  level.  By  mechanism  we  mean  the 
combination  of  conditions  that  result  in  a  tissue  level  accumulation  of  fracture.  In¬ 
trinsically,  the  micro  structure  of  bone  tissue  would  need  to  exceed  a  critical  strain 
state.  For  a  volume  that  is  deformed  quasi-statically,  the  critical  state  would  occur 
in  a  region  of  the  structure  that  experienced  the  largest  strains  due  to  geometric  or 
structural  features,  e.g.,  a  sharp  curve  or  point.  The  specific  location  of  the  fracture 
would  thus  largely  depend  on  the  structural  details  of  the  entire  bone  and  the  nature 
of  the  load  that  was  applied.  Alternatively,  the  material  may  fracture  due  to  a  com¬ 
bination  of  dynamic  effects,  such  as  elastic  waves  propagating  through  the  struc¬ 
ture,  inertial  loading,  or  the  extreme  case  of  a  shock  loading.  Ignoring  the  extreme 
case  of  shock,  any  dynamic  fracture,  or  more  accurately  large  strain  states,  within 
our  model  likely  originates  from  a  combination  of  inertial  and  wave  propagative  ef¬ 
fects.  Rate-dependent  constitutive  models,  whether  in  the  bone  or  in  the  flesh,  will 
change  transmissive  properties  in  the  frequency  domain;  however,  the  question  re¬ 
mains  whether  these  changes  are  large  enough  to  alter  the  observed  injury  location. 

4.3  Constitutive  Details  and  Loading 

We  return  to  2  assumptions  that  warrant  additional  discussion:  the  idealized  loading 
condition  and  the  use  of  isotropic  elastic  constitutive  models. 

The  loading  condition  applied  to  the  lower  leg  model  was  originally  developed  for 
Kolsky  bar  simulations  and  is  well  documented  in  Fitzpatrick  and  Scheidler.18  In 
their  work,  Fitzpatrick  and  Schiedler  concluded  that  piecewise  linear  functions  with 
discontinuities  in  the  derivative  of  applied  velocity  resulted  in  large  amounts  of 
ringing  and  high-frequency  noise.  Thus  a  high  quality  simulation  result  depended 
on  a  smooth  loading  function  with  continuity  in  the  accelerative  terms  as  well  as 
the  velocity.  Considering  the  discussion  of  layered  media  in  the  previous  section, 
and  as  we  move  toward  rate-dependent  constitutive  models,  control  over  the  high- 
frequency  content  and  noise  that  might  propagate  into  the  models  is  highly  desir¬ 
able.  We  note  that  the  largest  amplitude  high  frequency  content  of  our  loading  curve 
corresponds  to  frequencies  250  Hz  and  sub-kilohertz  (time  intervals  on  the  order  of 
milliseconds)  which  is  at  least  3  orders  of  magnitude  lower  than  the  periodic  behav¬ 
iors  of  the  transmission  coefficient  in  the  MHz  regime  discussed  in  Section  4.2. 
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The  results  from  Fig.  4  show  both  early  time  (less  than  1  ms)  and  late  time  (larger 
than  1  ms)  dependence  from  specific  elements.  Oscillations  are  observed  and  are 
most  prominent  in  the  bones  nearest  the  steel  plate  and  the  frequency  of  which  de¬ 
pended  on  the  material  parameters  of  the  flesh.  These  oscillations  at  early  times 
have  relatively  high  frequency  content  and  are  small  in  amplitude  when  compared 
to  the  lower  frequencies  that  are  larger  in  amplitude  at  later  times.  Although  fur¬ 
ther  analysis  is  needed,  we  can  speculate  that  the  changes  in  soft  tissue  properties 
quantitatively  change  the  transmission  coefficient  for  the  layer  of  cortical  bone  and 
selectively  suppress  portions  of  the  initial  transient  loading.  The  lower  frequencies 
at  later  times  are  essentially  seen  as  a  zero-frequency  component  and  are  largely 
transmitted  to  the  bone  regardless  of  material  parameter.  However,  if  we  are  miss¬ 
ing  a  key  feature  of  the  loading  by  assuming  an  idealized  function,  these  arguments 
would  no  longer  hold.  For  instance,  a  high-frequency  oscillation  superimposed  on 
our  loading  condition  might  transmit  extraordinarily  well  under  one  region  of  the 
calcaneus  compared  to  another  thicker  region.  Similarly,  these  arguments  may  not 
hold  at  large  strain  or  in  the  case  of  a  highly  viscoelastic  material.  For  instance,  a 
nonlinear  elastic  constitutive  model  for  flesh  would  predict  a  stiffened  elastic  mod¬ 
uli  in  addition  to  features  introduced  from  the  thickness  of  the  flesh  layer  changing. 
Another  question  that  arises  is  whether  the  loading  we  impose  is  overly  damaging. 
The  loading  is  entirely  specified,  which  potentially  delivers  an  unphysical  impulse. 

Typical  models  of  the  bone  have  included  elastic,24'29  elastic-plastic,1’21'22’30  and  vis¬ 
coelastic,20'31'32  along  with  various  failure  models  and  fracture  methods.  While  the 
loading  curves  may  replicate  the  rate-dependence  from  experiments,  the  unloading 
curves  are  rarely  taken  into  account.  For  a  plastic  material  this  return  is  quite  dif¬ 
ferent  than  a  viscoelastic  material.  Incorporating  the  wrong  physics  might  simply 
amount  to  introducing  the  same  error  as  not  including  the  rate-dependence  in  the 
first  place,  but  this  will  only  be  clear  through  detailed  comparisons  of  the  simula¬ 
tions  that  consider  both. 

There  is  utility  in  assuming  the  small  bones  of  the  lower  leg  are  undergoing  a  quasi¬ 
static  deformation.  This  is  supported  by  Figs.  2,  3,  and  4  along  with  the  estimates 
made  earlier  in  the  discussion.  If  this  is  the  case,  the  structures  of  these  bones  and 
the  nature  of  the  loading  would  be  most  critical  to  properly  understand  their  failure. 
Once  the  elastic  energy  enters  the  cortical  layer,  it  largely  follows  the  cortical  shell 
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since  waves  see  fairly  large  impedance  mismatches  between  the  flesh  and  cortical 
bone.  This  impedance  mismatch  can  delay  the  peak  of  loading,  similar  to  what  is 
observed  in  Fig.  4,  and  thus  play  a  role  in  the  distribution  of  stress  within  specific 
bones. 

In  past  and  current  US  Army  conflicts,  lower  extremity  injuries  amount  to  more 
than  50%  of  combat  wounds.7^9  The  prevalence  of  complications,  in  addition  to 
the  greater  disability  scores  seen  in  patients  with  foot  and  ankle  injuries,  stress  the 
importance  of  understanding  exactly  which  bones  are  injured.  If  a  rate-dependent 
model  can  change  the  location  of  where  we  see  injury,  we  need  to  understand  and 
incorporate  it  into  our  model.  The  dominant  features  we  need  to  understand  is  the 
inertial  loading,  the  wave  propagation  that  occurs  in  the  multilayered  structures, 
and  the  natural  structural  weaknesses  that  exist  in  each  of  the  bones.  The  latter  of 
which  can  be  understood  using  component  level  models. 

5.  Conclusion 

Injury  in  the  human  body  can  occur  over  multiple  length  and  time  scales.  It  links 
continuum  mechanics  to  biological  processes.  Simulations  and  other  modeling  ap¬ 
proaches  necessarily  simplify  the  problem  via  assumptions.  However,  significant 
research  questions  remain,  despite  (and  in  light  of)  these  simplifying  assumptions: 
What  physical  effects  carry  the  highest  order  terms  and  have  we  made  the  right  as¬ 
sumptions?  Does  injury  occur  on  a  strain  and  time  scale  where  nonlinearity  and 
rate-dependencies  are  large  effects  or  can  these  effects  be  neglected?  Where  should 
we  focus  our  attention;  to  the  variability  in  the  structures  and  interfaces  or  to  the 
constitutive  modeling  of  rate-dependent  and  nonlinear  materials?  In  conclusion, 
from  the  results  and  issues  brought  up  in  the  discussion,  it  is  imperative  that  we  re¬ 
visit  the  structural  and  interface  effects  within  our  model. 
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Appendix  A.  Details  of  the  Numerical  Simulation 
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Surface  geometries  for  the  anatomical  mesh  were  obtained  from  Zygote  Media 
Group,  Inc.  in  Wavefront  .obj  format.  These  surface  geometries  were  created  and 
“typified”  from  CT  (computer  tomography)  and  MRI  (magnetic  resonance  imag¬ 
ing)  scans  of  human  male  subjects  who  anthropomorphically  fit  within  the  50th  per¬ 
centile.  The  geometries  used  for  this  model  include  the  exterior  surfaces  of  the  skin, 
the  26  bones  of  the  foot,  the  patella,  and  the  long  bones  of  the  leg.  The  bones  and 
skin  were  initially  in  the  standing  position.  To  make  the  seated  position,  the  femur 
was  rotated  90°  in  the  sagittal  plane  and  then  visually  adjusted  in  the  frontal  and 
transverse  planes  so  that  the  femoral  condyles  sat  nicely  on  the  tibial  plateau.  The 
position  of  the  patella  was  also  adjusted  in  similar  fashion.  The  skin  required  man¬ 
ual  sculpting  and  patching  to  bend  at  the  knee  and  hips. 

Differentiation  between  cortical  and  trabecular  bone  types  was  achieved  by  dupli¬ 
cating,  downscaling,  manually  sculpting,  and  nesting  bone  geometries  to  create  the 
cortical  layer;  the  void  space  between  the  inner  scaled  bone  and  the  original  outer 
bone  was  considered  the  cortical  shell.  Thus,  a  sharp  interface  exists,  as  there  is  no 
transition  zone  from  cortical  to  trabecular  bone.  Appropriate  relative  cortical  thick¬ 
ness  was  sought  between  the  diaphysis  and  epiphyses  of  the  long  bones  (i.e.,  thicker 
in  the  diaphysis,  thinner  in  the  epiphyses).  However,  no  data  were  used  to  substan¬ 
tiate  this  effort.  Also,  the  inner  bone  structure  of  the  patella,  talus,  and  calcaneus  is 
not  well  understood,  and  therefore  a  uniformly  thick  cortical  shell  was  used.  The 
model  also  does  not  include  geometry  for  a  medullary  cavity;  the  interior  of  each 
bone  below  the  cortical  layer  in  its  entirety  is  considered  to  be  the  trabecular  bone 
region.  Bones  in  the  mid-and  forefoot  do  not  have  cortical-trabecular  bone  differ¬ 
entiation  and  are  assigned  cortical  bone  properties. 

All  manipulation  described  to  this  point  was  done  in  Blender  (V2.49,  Amster¬ 
dam,  Netherlands).  Stereolithography  (STL)  files  were  exported  from  Blender  and 
brought  into  the  meshing  software  ICEM  CFD  (V15;  ANSYS,  Canonsburg,  PA). 
Further  geometry  manipulation  was  performed  in  ICEM  CFD;  the  middle  and  dis¬ 
tal  phalanges  of  the  toes  were  shrink-wrapped  so  they  were  fused  together.  Every¬ 
thing  was  then  meshed  simultaneously,  hence  the  nodes  between  different  compo¬ 
nent  parts  are  shared  (i.e.,  between  cortical  and  trabecular  bone,  cortical  bone,  and 
surrounding  flesh). 
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The  model  was  then  exported  from  ICEM  CFD  in  LS-Dyna  format  and  opened  in 
LS-PrePost  (V4.0,  LSTC,  Lawrence  Livermore  National  Laboratory).  Mesh  manip¬ 
ulation  was  performed  here,  which  included  cleaning  up  stray  elements  that  would 
connect  neighboring  bones;  these  elements  were  manually  selected  and  reassigned 
to  the  surrounding  flesh.  Hence,  each  bone  in  the  model  is  isolated  from  any  sur¬ 
rounding  bones  via  a  flesh  buffer.  A  flat  plate  of  hexahedral  elements  was  created 
in  LS-PrePost  and  centered  under  the  feet.  The  model  was  then  initialized  in  LS- 
Dyna  version  6.1.0  to  bring  the  bottom  of  the  feet  flat  and  in  contact  with  the  plate. 
This  included  applying  a  velocity  of  0.5  m/s  to  the  plate  in  the  upward  (Z)  direction 
toward  the  feet.  Once  the  bottom  of  the  heel  contacted  the  plate,  the  nodal  coordi¬ 
nates  of  the  model  at  that  state  were  saved  and  exported.  Next,  a  new  instance  of  the 
original  model  was  opened  in  LS-PrePost,  and  the  initialized  state  coordinates  were 
imported  without  offset  to  rewrite  the  nodal  coordinates  of  the  original  model.  For 
this  work,  the  left  leg  and  pelvis  were  removed  from  the  model,  leaving  the  right 
leg.  The  flesh  was  removed  just  above  midfemur;  the  entire  right  femur  itself  re¬ 
mained  intact.  A  jagged  mesh  at  the  site  of  flesh  deletion  remained  and  was  cleaned 
up  by  hand  by  creating  elements  from  nodes  on  the  exposed  face  of  the  flesh. 

To  run  the  model  in  Sierra/S olidMechanics,  the  mesh  needed  to  be  in  Exodus  file 
format,  which  LS-PrePost  cannot  directly  export.  Hence,  the  model  was  converted 
to  an  Exodus  file  format  through  the  use  of  HyperMesh  (VI 2.0;  Altair)  and  Cubit 
(V13.1;  Sandia  National  Laboratory);  the  model  was  exported  in  LS-Dyna  format 
from  LS-PrePost,  imported  into  HyperMesh,  exported  as  an  Abaqus  .inp  file,  im¬ 
ported  into  Cubit,  and  then  exported  as  an  Exodus  .g  file.  Tied  contact  was  used  be¬ 
tween  the  flesh  and  bones. 

The  model  contains  approximately  a  half  million  tetrahedral  elements.  The  typi¬ 
cal  criticism  of  the  tetrahedron  element  is  volumetric  locking  for  materials  with  a 
large  Poisson’s  ratio.  The  large  values  of  Poisson’s  ratio  of  the  flesh  would  cer¬ 
tainly  cause  locking  in  the  traditional  tet  element.  As  a  result,  we  are  using  a  nodal- 
based  tet1.  For  bend  tests  this  nodal-based  tet  did  not  suffer  from  volumetric  lock¬ 
ing  for  Poisson’s  ratios  as  high  as  0.49.  While  they  do  not  always  perform  as  well 

'Dohrmann  C,  Heinstein  M,  Jung  J,  Key  S,  Witkowski  W.  Node -based  uniform  strain  ele¬ 
ments  for  three-node  triangular  and  four-node  tetrahedral  meshes.  International  Journal  for  Numer¬ 
ical  Methods  in  Engineering.  2000;47(9):  1549-1568. 
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as  a  traditional  hexahedral  element,  given  the  complicated  geometry  of  the  flesh, 
they  were  the  best  available  option.  A  rate  integrated  formulation  of  a  linear  elastic 
material  was  used  in  SIERRA  for  all  anatomical  components.  Typical  stresses  ob¬ 
served  in  the  bones  were  small  compared  to  the  elastic  moduli,  thus  implying  that 
strains  were  also  small. 

Postprocessing  of  simulation  results  was  carried  out  in  Para  View  (V3.14.0;  Kit- 
ware),  SEACAS  (Sandia  National  Laboratory),  EnSight  (vl0.0.3c;  CEI,  Inc.,  Apex, 
NC),  and  MATLAB  (The  MathWorks,  Natick,  MA). 
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Appendix  B.  Mathematical  Formulation  of  a  Smooth  Loading  Pulse 
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For  all  of  the  simulations  of  the  lower  leg,  we  require  a  loading  pulse  (i.e.,  an  applied 
velocity  history  v)  that  ramps  up  smoothly  to  a  constant  value  before  subsequently 
unloading  (Fig.  lc  in  report).  In  this  appendix  we  formulate  a  parameterized  loading 
pulse  using  a  Hermitian  smoothing  spline. 

We  use  a  Hermitian  smoothing  spline  of  degree  5  as  the  basis  for  the  ramp-up 
portion  of  the  loading  pulse1): 

H5(x)  =  10a;3  —  15a;4  +  6a;5,  0  <  x  <  1 .  (B-l) 


This  function  increases  from  0  at  x  —  0  to  1  at  x  —  1,  and  has  zero  first  and 
second  derivatives  at  both  of  these  end  points.  Following  Fitzpatrick  and  Scheidler1 
(without  loss  of  generality  letting  strain-rate  become  velocity  v  and  strain  becomes 
displacement  5),  one  can  use  H5(x)  to  construct  a  velocity  function  that  ramps  up 
from  0  to  Umax  over  the  time  interval  [  0  ,  tr  ],  maintains  this  constant  value  over  the 
time  interval  [tr  ,tr  +  tp],  and  then  unloads  to  zero  over  the  time  interval  [tr  + 


tp  , td  ]  • 


0, 


^max  H5  ( 


trn 


t  <  0 


0  <  t  <  tr 


v(t)  =  < 


^max: 

^max 

0, 


^  t—(tr+tp)  ^ 


t  ty  |  tp 
,  ty  tp  ^  t  ''' 

t>td. 


(B-2) 


tr  is  referred  to  as  the  rise  time  (i.e.,  the  time  it  takes  to  reach  a  constant  velocity), 
tp  is  the  duration  of  the  constant  velocity  plateau,  and  td  is  the  total  duration  of  the 
pulse.  In  the  simulations,  we  made  the  additional  assumption  that  the  loading  pulse 
is  symmetric;  therefore,  the  unloading  duration  td  —  (tr  +  tp)  equals  the  rise  time, 

so  that  td  =  2  tr  +  tp. 

'Fitzpatrick  J,  Scheidler  M.  Considerations  for  numerical  modeling  of  dynamic  behavior  of  soft 
materials.  In:  Chalivendra  V,  Song  B,  Casern  D,  editors.  Dynamic  Behavior  of  Materials,  Volume  1. 
Proceedings  of  the  2012  Annual  Conference  on  Experimental  and  Applied  Mechanics;  2012  June 
11-14;  Costa  Mesa,  CA.  New  York  (NY):  Springer;  2012. 
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Note  that 


'y(O)  =  v(tr)  =  v(tr  +  tp)  =  0 


(B-3) 


and 


b(0)  =  v(tr)  =  v(tr  +  tp)  —  0  . 


(B-4) 


Since  v  and  v  are  identically  zero  during  the  constant  velocity  plateau,  as  well  as 
for  t  <  0  and  t  x  / ( / .  the  conditions  in  Eqs.  B-3  and  B-4  guarantee  that  v  is  twice 
continuously  differentiable.  This  degree  of  smoothness  was  necessary  to  avoid  an 
unwanted  ringing  in  the  simulations  that  would  otherwise  result  from  a  lower  order 
polynomial. 


An  explicit  expression  for  S(t )  is  rather  messy,  and  in  practice  it  is  simpler  to  work 
with  the  integral  of  H5: 

H5(x)  =  j_‘  H5{y)  d  y  =  x4  (x2  -  3*  +  0  .  (B-5) 

Integrating  each  case  in  Eq.  B-2  with  respect  to  time,  assuming  a  symmetric  loading 
pulse,  and  using  IHI5 (1)  =  1/2,  gives: 


m 


0, 

^max  tr  H5  (  ^  )  ? 

^max  [tr / 2  +  (t  —  tr)] , 
tr/2  +  (t  -  tr)  - 

^max  [tr  +  tp] , 


t  <  0 
0  <  t  <  tr 

tp  ^  t  <\  tp  |  tp 

tr  “t-  tp  ^  t  fC  2  tr  +  tp 

t  ~^>  2 \tp  |  tp  —  t d  • 

(B-6) 
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Intentionally  left  blank. 
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Appendix  C.  Failure  Criteria  Figures 
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Failure  occurrences  are  tallied  at  various  bones  for  a  range  of  thresholds.  Each 
panel  is  specific  to  an  individual  bone  and  tallies  the  number  of  simulations  for 
which  the  absolute  value  of  the  pressure  within  any  element  in  the  bone  exceeded 
the  damage  threshold.  The  bar  colors  indicate  which  damage  threshold  the  tally 
corresponds  to.  Low  and  high  values  for  the  material  parameter  as  well  as  the  spe¬ 
cific  damage  thresholds  used  are  indicated  on  the  right  under  the  key.  Figures  C- 
1  through  C-5  look  at  a  pressure-based  failure  criteria.  Figures  C-6  through  C-10 
look  at  a  von  Mises-based  failure  criteria. 


Calcaneus  Talus  Tibia  Fibula  Forefoot  Key 


Fig.  C-l  Excessive  pressure  failure  occurrence  seen  in  each  bone  by  varying  the  Young’s  mod¬ 
ulus  of  cortical  bone 
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Fig.  C-2  Excessive  pressure  failure  occurrence  seen  in  each  bone  by  varying  the  density  of 
trabecular  bone 
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Fig.  C-3  Excessive  pressure  failure  occurrence  seen  in  each  bone  by  varying  the  Young’s  mod¬ 
ulus  of  trabecular  bone 
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Fig.  C-4  Excessive  pressure  failure  occurrence  seen  in  each  bone  by  varying  the  Poisson’s  ratio 
of  trabecular  bone 
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Fig.  C-5  Excessive  pressure  failure  occurrence  seen  in  each  bone  by  varying  the  Poisson’s  ratio 
of  flesh 
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Fig.  C-6  Excessive  von  Mises  failure  occurrence  seen  in  each  bone  by  varying  the  Young’s 
modulus  of  cortical  bone 


Calcaneus  Talus  Tibia  Fibula  Forefoot  Key 


Fig.  C-7  Excessive  von  Mises  failure  occurrence  seen  in  each  bone  by  varying  the  density  of 
trabecular  bone 
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Fig.  C-8  Excessive  von  Mises  failure  occurrence  seen  in  each  bone  by  varying  the  Young’s 
modulus  of  trabecular  bone 
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Fig.  C-9  Excessive  von  Mises  failure  occurrence  seen  in  each  bone  by  varying  the  Poisson’s 
ratio  of  trabecular  bone 
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Fig.  C-10  Excessive  von  Mises  failure  occurrence  seen  in  each  bone  by  varying  the  Poisson’s 
ratio  of  flesh 
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